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Abstract 

The  thrust  of  this  project  has  been  two-fold,  namely,  (1)  to  improve  the 
performance  of  existing  models  of  Earth’s  gravitational  field,  mainly  with 
respect  to  speed  of  evaluation;  and  (2)  to  develop  new  multiresolution  esti¬ 
mation  techniques  to  produce  new  gravitational  models. 

We  have  constructed  two  local  models  of  the  gravitational  field,  both  of 
which  map  the  surface  of  a  sphere  to  the  surface  of  a  cube.  These  models 
differ  in  the  choice  of  basis  functions.  The  first  uses  multiwavelets  to  represent 
the  gravitational  field  at  a  fixed  distance  from  Earth,  and  the  second  model 
uses  B-splines.  Both  models  use  polynomial  interpolation  to  compute  the 
variation  in  height  of  the  gravity  field. 

Significant  progress  has  been  made  towards  implementing  a  procedure 
for  estimating  gravity  models  directly  from  physical  measurements.  We  have 
developed  a  new  estimation  algorithm,  a  multiresolution  rank-revealing  QR 
decomposition.  This  algorithm  produces  the  “minimum  detail”  solution  in¬ 
stead  of  the  usual  minimum  norm  solution  of  the  ill-conditioned  least  squares 
problem. 

A  basis  for  bandlimited  functions  has  been  constructed  using  a  new 
method  for  computing  the  generalized  Gaussian  quadratures  for  exponen¬ 
tials.  These  bases  are-closely  related  to  the  prolate  spheroidal  wave  functions, 
and  we  plan  to  create  the  next  generation  of  models  for  evaluation  and  esti¬ 
mation  of  the  gravity  field  using  such  bases.  These  bases  are  nearly  optimal 
in  terms  of  the  number  of  coefficients  necessary  to  represent  a  bandlimited 
function. 

In  addition,  we  have  begun  work  on  a  new  type  of  ODE  solver  (spectral 
deferred  corrections),  which  will  work  in  conjunction  with  the  new  gravity 
models,  and  should  provide  further  improvements  to  speed  and  accuracy  of 
computing  satellite  ephemerides. 

New  gravity  models  of  the  type  described  here  were  transferred  to  Space 
Warfare  Command  in  Colorado  Springs. 


1  Introduction 


We  recall  the  form  of  the  spherical  harmonic  model  of  the  gravitational  po¬ 
tential 

nr,(6,»)  =  ^|i  +  f;(f)"ms)| ,  a) 

where  GM  is  Earth’s  gravitational  constant,  r  is  the  length  of  the  radius 
vector  from  Earth’s  center  of  mass,  R  is  the  radius  of  the  Earth  (either 
equatorial  or  polar  radius),  4>  is  geocentric  longitude  and  9  is  geocentric 
latitude.  The  spherical  harmonic  Yn((j),9)  is  defined  as 

Yn{<f>,  9)  =  (sin  9)  ( C ™  cos  m<f>  +  sin  m</>)  , 

0 

where  C°n,  C\, . . . ,  C^,  5* , . . . ,  S%,  for  2  <  n  <  N  are  normalized  coefficients, 
and  is  the  normalized  associated  Legendre  function  of  degree  n  and  order 
m.  The  parameter  N  in  (1),  the  number  of  cerms  retained  in  the  model, 
determines  the  order  of  the  model.  For  example,  N  =  41  corresponds  to  a 
41st  order  model,  and  so  on.  As  is  well  known,  (1)  is  a  solution  of  the  Laplace 
equation  in  spherical  coordinates  (r,4>,9),  r  >  R 

The  cost  of  evaluating  V  at  a  point  (r,  <p,  9)  via  (1)  grows  rapidly  with  the 
number  of  retained  terms.  Namely,  if  N  terms  are  retained  in  (1)  then  the 
number  of  operations  to  evaluate  V  is  proportional  to  N2.  Thus,  changing 
the  model  from  that  using  20  terms  to  the  model  using  200  terms  requires 
roughly  100  times  more  operations. 

We  observe  that  there  are  several  sources  of  inefficiency  in  using  this  rep¬ 
resentation.  First,  the  functions  ( R/r)n ,  n  =  2,3,...  are  “nearly  parallel,” 
which  means  that  the-organization  of  the  sum  in  (1)  is  not  efficient  for  com¬ 
puting.  Second,  the  functions  Yn(<t>,9)  are  global  and  thus  cannot  model  the 
regional  variations  in  the  Earth’s  geopotential  with  the  detail  and  economy 
provided  by  local  functions  such  as  splines  or  wavelets.  Unlike  global  func¬ 
tions,  the  resolution  of  local  functions  can  be  adjusted  to  the  level  required 
to  accurately  model  local  features,  a  process  known  as  “adaptive  sampling.” 

In  what  follows  we  use  the  terms  “local”  and  “multiresolution”  representa¬ 
tions.  In  a  local  representation  the  elementary  building  blocks  are  basis  func¬ 
tions  with  localized  support,  e.g.  B-splines.  These  representations  typically 
do  not  use  adaptive  sampling;  rather,  the  resolution  is  uniform  throughout 
the  model.  Multiresolution  representations  also  use  localized  basis  functions 
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(e.g.  wavelets)  but,  in  addition,  can  use  adaptive  sampling  to  more  efficiently 
accommodate  different  variability  of  the  geopotential  in  different  regions. 


2  Local  Models  of  the  Gravity  Field 

Two  local  models  have  been  developed  under  this  grant,  but  before  describing 
them  it  will  be  advantageous  to  provide  background  by  describing  an  earlier 
model,  developed  under  an  earlier  grant. 

2.1  Doubly  Periodic  B-spline  Model 

The  local  model  developed  under  DARPA  grant  Efficient  Representation 
of  the  Earth’s  Gravitational  Field  comprises  a  set  of  spherical  shells,  each 
corresponding  to  a  fixed  value  of  the  radial  distance  r.  On  each  shell,  the 
potential  field  (or  one  of  its  derivatives)  is  modeled  by  a  doubly  periodic 
B-spline  expansion.  For  purposes  of  illustration,  let  ro  be  a  fixed  value  of  r. 
To  model  the  field  on  a  spherical  shell  at  r  =  r0  we  construct  the  following 
representation, 


M—l  M-l 

V(r0,<f>,9)  =  Sk,tP(Mx-k)p(My-l)  (2) 

*=o  1=0 

where  2-kx  =  9  and  2iry  =  <f>.  In  (2),  /3  denotes  the  B-spline  (compactly 
supported,  piecewir„  polynomial  function)  of  sufficiently  high  order,  and  1/M 
is  the  length  of  th ?.  largest  interval  on  which  the  spline  is  a  polynomial.  The 
choice  of  the  integer  M  and  the  order  of  the  B-splines  depends  upon  accuracy 
and  memory  requirements  and  both  are  adjustable  parameters  in  the  model. 

The  cost  of  evaluating  the  B-spline  series  (2)  is  a  constant  that  depends 
on  the  order  of  the  spline,  and  is  proportional  to  ( mdeg  +  l)2,  where  mdeg 
denotes  the  degree  of  the  spline.  Thus,  for  example,  if  we  use  B-splines  of 
degree  seven,  then  roughly  sixty-four  multiplications  are  required  to  evaluate 
the  right-hand  side  of  (2).  Note  that  this  cost  is  independent  of  the  parameter 
M ,  which  instead  governs  memory  storage  requirements,  with  its  magnitude 
being  dictated  by  the  required  precision.  It  is  obvious  from  (2)  that  M2 
coefficients  must  be  stored  in  computer  memory  for  each  spherical  shell.  The 
values  of  M  and  mdeg  are  in  inverse  proportion;  the  larger  is  M,  the  smaller 
is  mdeg ,  and  vice  versa,  in  order  to  achieve  any  prescribed  accuracy.  Thus, 
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the  issue  of  parameter  selection  becomes  a  trade-off  between  computational 
speed  and  computer  storage  requirements. 

Now  let  us  describe  how  we  use  (2)  to  obtain  the  values  of  the  geopotential 
and  its  derivatives  at  locations  between  the  tabulated  shells.  Suppose  that 
we  want  to  evaluate  at  a  distance  r'  above  the  Earth’s  surface,  and  that 
rk  <  r'  <  rk+ 1,  where  rk  and  rk+\  correspond  to  consecutive,  tabulated 
shells.  We  simply  compute  the  value  of  the  potential  on  the  shells  at  rk 
and  rk+i  at  the  proper  locations  in  ( <t>,9 ),  then  construct  the  interpolating 
polynomial  that  joins  these  two  points  and  evaluate  this  polynomial  at  r  =  r'. 
This  is  the  procedure  for  a  linear  interpolating  polynomial.  In  practice,  we 
ordinarily  use  a  higher  order  polynomial,  for  example  fifth  degree,  which 
requires  evaluation  on  more  spherical  shells  but  the  procedure  is  otherwise 
identical. 

The  total  cost  of  evaluation  per  point  for  the  B-spline  model  can  now  be 
made  precise.  If  the  degree  of  the  B-spline  used  is  mdeg,  and  the  degree  of 
the  polynomial  used  to  interpolate  between  shells  is  ideg,  then  the  total  cost 
of  evaluating  the  model  at  a  point  is  C  *  ( ideg  +  1)  *  ( mdeg  +  l)2,  where  C 
is  a  constant.  Note  that  ideg  is  also  an  adjustable  parameter  in  this  model, 
and  its  value  can  be  made  smaller  if  we  are  willing  to  make  M  larger,  as 
described  above  for  the  parameter  mdeg. 

We  now  outline  the  steps  in  computing  the  coefficients  for  the  B-spline 
expansion  (2).  For  efficient  and  accurate  computation  of  the  coefficients,  it 
is  convenient  to  extend  the  potential  function  V  so  that  it  is  27r-periodic  in 
both  9  and  <f>.  This  is  done  by  defining 

f  V (r,  <f>  +  7T,  —7 r  —  9) ,  for  — 7r  <  6  <  — 7r/2 
Vp(r,  ct>,6)  =  \  V(r,  <t>,  9)  ,  for  -tt/2  <9  <  tt/2 

[  V(r, <f>  +  7r, 7r  —  9)  ,  for  7r/2  <  0  <  7r. 

for  9  e  [ — 7r, 7r]  ( V  is  already  27r-periodic  in  0). 

There  are  several  useful  bases  for  splines,  in  addition  to  the  “canonical” 
B-splines,  and  in  this  model  generation  step  we  also  use  the  interpolating 
splines,  L(x).  Interpolating  splines  satisfy  L(0)  =  1,  and  L(n )  =  0  if  n  is  a 
non-zero  integer.  This  property  enables  us  to  first  write  the  spline  expansion 


as 


V(r0,<t>,9) 


£  AwL 


k=0  1=0 


with  Akii  =  V(r0,<j)i,9k),  where  9k  =  2-irk/M  and  <f>t  =  2i d/M.  To  change 
the  basis  from  interpolating  splines  (which  are  difficult  to  evaluate)  to  the 
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B-splines,  all  we  have  to  do  is  apply  an  FFT  to  the  coefficient  matrix 
modify  the  Fourier  transform  by  a  factor,  then  apply  the  inverse  FFT  to 
obtain  the  B-spline  coefficients  in  (2).  As  a  result,  the  representation  (2) 
interpolates  the  potential  function  V(r.  <jr,  0),  for  r  =  r0,  on  the  equispaced 
grid  (2irk/M,  2i xl/M). 

Once  the  representation  (2)  for  the  geopotential  is  obtained,  it  is  easy 
to  build  its  reduced  version,  which  is  obtained  by  simply  removing  high 
frequencies  from  the  model  by  application  of  multiresolution  decomposition 
to  the  B-spline  expansions.  To  describe  the  construction  process,  let  us 
assume  that  M  =  2-J,  where  j  <  0  is  an  integer,  and  denote  the  coefficients 
in  (2)  by  (4  J,  to  indicate  that  they  belong  to  the  scale  j. 

The  result  of  applying  one  step  of  multiresolution  decomposition  is  to 
compute  the  coefficients  {4]/1}  on  the  next  coarser  scale  (j  +  1)  (i.  e.  the 
scale  corresponding  to  double  the  original  step  size)  from  the  coefficients  on 
scale  j.  From  M  coefficients  per  row  on  scale  j  we  obtain  M/2  =  2“t?+1) 
coefficients  per  row  on  scale  ( j  +  1).  The  steps  are  as  follows: 

1.  Apply  an  FFT  to  the  coefficient  matrix  {4/}  to  obtain  the  matrix 

RJ- 

2.  Apply  a  one-dimensional  decomposition  in  each  index,  i.  e.  first  on  rows, 
then  on  columns,  of  the  matrix  (4,J-  The  one-dimensional  decompo¬ 
sition  is  defined  as  follows, 

=  \  {*0  (^)  Sm  +  *0  +  it)  4+w} 

for  m  =  0, . . . ,  M/2  —  1,  where  the  function  fho  is  a  27r-periodic  function 
associated  to  the  B-spline. 

3.  Apply  an  inverse  FFT  to  the  matrix  {4^}  to  obtain  the  coefficient 
matrix  {4^}  on  the  coarser  scale. 

Removal  of  the  higher  frequencies  from  the  model  in  this  manner  is  more 
stable  than  simply  truncating  them  as  is  done,  for  example,  in  truncation  of 
the  WGS84-70  model  to  obtain  the  WGS84-41  model  [11].  If  desired,  Steps 
1-3  above  can  be  repeated,  to  obtain  coefficients  (4^2}>  so  on. 
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2.2  Multiwavelet  Cube  Model 

In  this  model,  the  surface  of  the  sphere  is  mapped  to  the  surface  of  a  cube 
and,  instead  of  spheres,  the  concentric  shells  form  a  sequence  of  nested  cubes. 
A  point  on  the  surface  of  a  sphere  is  mapped  to  a  point  on  the  “reference 
cube”  (which  has  faces  perpendicular  to  the  coordinate  axes  and  at  a  distance 
of  one  unit  from  the  origin)  using  the  following  simple  algorithm. 

•  Input  coordinates  (r,  <f>, 0)  on  the  spherical  surface  of  radius  r. 

•  Compute  (x  =  r sin# cos 0,  y  =  r sin 9 sin <j>,  z  =  rcos9). 

•  Find  d  =  max{|x|,  |y|,  |z|}. 

•  Coordinates  on  the  reference  cube  are  (£,77,  C)  =  (x/ d,  y/ d,  zj d). 

Geometrically,  we  can  think  of  a  ray  emanating  from  the  origin,  and  inter¬ 
secting  the  sphere  and  the  reference  cube  each  at  a  single  point.  These  two 
points  are  then  mapped  one  to  the  other. 

We  next  place  an  equispaced  rectangular  grid  on  each  face  of  the  cube, 
which  can  then  be  mapped  backwards  onto  the  sphere.  We  note  that  dis¬ 
tortion  of  the  grid  on  the  spherical  surface  caused  by  the  curvature  of  the 
sphere  is  limited,  and  does  not  cause  any  problem  near  the  poles. 

The  rectangular  grid  partitions  each  face  of  the  cube  into  a  number  of 
square  subdivisions,  and  we  build  a  wavelet  representation  of  the  geopoten¬ 
tial  on  each  subdivision.  Currently,  the  basis  functions  we  use  are  the  mul¬ 
tiwavelets  (see  [1]  and  [2]),  chosen  because  they  form  an  orthonormal  basis 
on  a  square  subdivision  without  overlapping  into  adjacent  subdivisions.  . 

The  scaling  functions  in  a  multiwavelet  basis  are  classes  of  orthogonal 
polynomials,  for  example  Legendre  polynomials.  Formulas  for  interpolat¬ 
ing  the  geopotential  at  the  Gaussian  nodes  within  the  subinterval  are  easily 
obtained  using  the  Gaussian  integration  rules. 

To  evaluate  the  field  between  tabulated  shells,  we  construct  the  La- 
grangian  polynomial  that  interpolates  the  values  on  several  adjacent  shells. 
Alternatively,  we  can  also  use  expansions  of  multiwavelets  to  represent  the 
variation  in  height  of  the  gravity  field. 


2.3  Spline  Cube  Model 

Concerning  performance  of  the  two  models  described  thus  far,  two  observa¬ 
tions  have  been  made: 
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1.  the  multiwavelet  cube  model  has  a  more  efficient  memory  access  than 
the  doubly  periodic  spline  model,  resulting  in  better  speed  for  evalua¬ 
tion; 

2.  the  doubly  periodic  spline  model  requires  substantially  fewer  coeffi¬ 
cients  to  achieve  the  same  accuracy  as  the  multiwavelet  cube  model. 

The  model  described  in  this  section  has  been  successful  in  combining  the  best 
features  of  both  these  models. 

To  evaluate  our  models  it  is  necessary  to  interpolate  between  the  tabu¬ 
lated  shells.  To  do  this  we  must  evaluate  the  two-dimensional  expansions 
at  a  given  point  on  several  consecutive  shells.  However,  the  computation 
involves  only  a  few  coefficients  on  each  shell,  and  speed  is  gained  in  memory 
access  if  these  can  be  stored  closer  together. 

This  can  be  done  by  subdividing  the  surface  into  a  number  of  “panels” 
which,  taken  together,  cover  the  spherical  surface.  Panels  for  the  same  angles 
but  different  heights  are  stored  contiguously  in  memory,  which  allows  for 
faster  memory  access. 

In  this  model,  a  partition  of  the  spherical  surface  is  accomplished  by 
subdividing  each  spherical  shell  into  six  panels,  which  may  be  regarded  as 
the  six  faces  of  a  cube.  This  is  done  in  such  a  way  that  the  grid  spacing 
for  the  polar  regions  is  the  same  as  that  for  the  equatorial  regions,  and  no 
excessive  distortion  of  the  grid  on  the  spherical  surface  occurs. 

To  obtain  six  square  panels,  we  subdivide  the  surface  of  a  sphere  as 
indicated  in  the  following  table. 


□ 

angle  ranges 

x-coordinate 

y-coordinate 

m 

—  7T  <(j)<  — 7r/2  , 

-it/ 4  <  #  <  tt/4 

a  (f)  +  3 

a# 

B 

— 7r/2  <  (j)  <  0  , 

— 7t/4  <  #  <  7t/4 

a  <f>  +  1 

a6 

a 

0  <(j><  7t/2  , 

— IT / 4  <  #  <  tt/4 

a  (j>  —  1 

a  9 

B 

7t/2  <  <t>  <  7T  , 

— 7t/4  <  #  <  7t/4 

acf>  —  3 

a# 

B 

M  <  1 »  M 

—a  sin-1  (7) 

K3 

ItI  <  1 ,  ^ 

>  1,  #  <  0 

—a  sin-1  (7) 

where  a  =  4/7T,  u>  —  tan#/ cos <f>,  and  7  =  cos# sin <j>.  Coordinates  on  the 
face  of  each  panel  are  (x,  y),  where  —1  <  x,y  <  1.  The  panels  designated  5 
and  6  contain  the  north  and  south  poles,  respectively.  For  the  x-coordinate 
in  panels  5  and  6,  we  use  the  minus  sign  if  u  >  0  and  the  plus  sign  if  a;  <  0. 
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We  note  that  the  B-spline  expansion  for  each  panel  overlaps  its  immediate 
neighbors.  Thus,  to  use  this  model  for  the  estimation  problem,  we  would  need 
to  add  a  certain  number  of  equations  to  ensure  that  the  representation  near 
the  boundaries  of  each  panel  matches  with  that  of  its  neighbors. 

Coefficients  for  this  model  are  computed  in  the  same  manner  as  described 
above  for  the  doubly  periodic  spline  model.  After  computing  coefficients 
for  a  doubly  periodic  B-spline  expansion  to  cover  the  sphere,  coefficients 
sufficient  to  cover  each  of  the  first  four  panels  are  extracted  and  stored.  The 
sphere  is  then  rotated,  another  doubly  periodic  expansion  is  computed,  then 
coefficients  are  extracted  to  cover  the  two  remaining  panels.  The  rotation  is 
necessary  to  avoid  excessive  distortion  of  the  grid  near  the  poles. 

We  note  that  the  gravity  force  field  for  the  polar  regions  must  be  stored 
in  rectangular  coordinates.  In  the  spherical  coordinate  system,  the  represen¬ 
tation  of  a  vector  at  one  of  the  poles  is  not  unique: .  the  components  of  the 
vector  can  vary  depending  on  the  orientation  of  the  spherical  unit  vectors. 
This  orientation  depends  on  the  angle  (f>,  which  is  indeterminate  at  the  poles. 
Thus,  although  the  vector  at  the  pole  is  unique,  its  projection  onto  spherical 
unit  vectors  is  not. 

2.4  Performance  Results 

The  table  below  contains  performance  results  related  to  the  doubly  periodic 
spline  model  and  the  spline  cube  model.  The  parameter  “model  order”  in 
the  far  left  column  refers  to  the  number  of  terms  retained  in  the  spherical 
harmonic  model.  The  next  column  contains  the  size  of  the  spherical  harmonic 
model  in  megabytes.  As  can  be  seen,  the  spherical  harmonic  model  requires 
very  little  storage  space.  The  next  pair  of  columns  contains  information 
relating  to  the  doubly  periodic  spline  model,  namely  the  size  in  megabytes 
and  speed-up  factor.  The  speed-up  factor  is  obtained  by  measuring  the  time 
required  to  compute  the  gravity  vector  at  a  large  number  of  points,  using 
both  the  spline  model  and  the  spherical  harmonic  model,  then  dividing  the 
time  for  the  spherical  harmonic  model  by  the  time  for  the  spline  model.  The 
final  pair  of  columns  contains  size  in  megabytes  and  speed-up  factor  for  the 
spline  cube  model.  (It  turns  out  that  the  multiwavelet  cube  model,  as  it  is 
currently  implemented,  requires  excessive  memory  storage,  and  we  chose  not 
to  report  performance  results  for  this  model  at  this  time.) 
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Model 

Order 

Size  of 
Sph  Model 

Doubly  Periodic 

Spline  Cube 

Size 

Size 

Speed-up 

18 

0.003 

4.65 

2.3 

10.3 

3.6 

41 

0.013 

33.10 

9.0 

27.9 

19.5 

70 

0.038 

76.48 

22.9 

70.3 

56.1 

The  WGS84  spherical  harmonic  model  [11]  was  used  for  these  tests,  with 
agreement  to  about  10~u  between  spherical  harmonic  and  spline  models 
being  maintained  throughout. 

Note  that  storage  requirements  for  both  spline  models,  while  much  larger 
than  what  is  needed  for  the  spherical  harmonic  models,  are  nevertheless 
quite  reasonable.  Furthermore,  the  spline  models  have  achieved  a  significant 
speed-up  in  evaluation  time  over  the  spherical  harmonic  models.  Note  also 
that  the  spline  cube  model,  which  was  developed  during  the  past  year,  has 
gained  better  than  a  factor  of  two  in  speed  over  the  previous  spline  model, 
in  addition  to  reducing  memory  requirements  somewhat. 


3  Estimation  of  the  Geopotential 

Let  us  explain  briefly  why  it  is  necessary  to  develop  new  estimation  al¬ 
gorithms.  The  process  of  estimation  (model  building)  of  the  geopotential 
presents  several  problems.  First,  the  spatial  frequency  bandwidth  is  chang¬ 
ing  both  as  a  function  of  the  distance  from  the  Earth  and  the  location  on 
the  Earth’s  surface.  Second,  the  data  are  typically  collected  on  an  unequally 
spaced  grid  and,  as  a  result,  the  condition  number  of  matrices  involved  in 
solving  the  estimation  problem  is  usually  very  high. 

Within  the  spherical  harmonic  representation,  the  problem  of  increasing 
the  accuracy  of  the  model  was  addressed  by  increasing  the  number  of  terms  in 
the  expansion.  The  difficulty  with  this  approach  is  that  spherical  harmonics 
(being  global,  oscillatory  functions)  depend  on  cancellation  (destructive  in¬ 
terference)  to  achieve  the  approximation — changing  even  a  single  coefficient 
in  the  model  has  a  global  effect.  In  consequence,  algorithms  for  model  gener¬ 
ation  are  expensive,  since  we  have  to  solve  systems  with  dense  (full)  matrices. 
It  is  difficult,  if  not  impossible,  to  adjust  the  spatial  frequency  contents  lo¬ 
cally.  In  particular,  there  is  a  real  problem  in  incorporating  observations  of 
the  gravitational  potential  near  the  surface  (e.g  on  the  ground)  with  those 
obtained  on  satellites,  apparently  due  to  the  different  spectral  contents  of 
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the  data. 

Our  goal  has  been  to  develop  multiresolution  models  of  the  Earth's  grav¬ 
itational  potential  that  do  not  suffer  from  the  difficulties  for  use  and  es¬ 
timation  outlined  above.  Multiresolution  models  use  basis  functions  (e.g. 
wavelets)  with  localized  support  in  both  space  and  spectral  domains.  This 
allows  us  to  generate  models  Where  changes  in  most  of  the  parameters  will 
produce  only  local  changes  in  the  geopotential  (up  to  any  finite  but  arbitrary 
accuracy) . 

For  a  gravitational  field  this  statement  might  sound  strange  since  such 
fields  always  give  rise  to  long  range  interactions;  in  particular,  there  are 
no  negative  masses  and  interaction  appears  to  be  global.  It  is  easier  to 
understand  this  point  using  the  analogy  with  electrostatic  fields  which  are 
mathematically  equivalent  but  allow  us  to  consider  negative  charges.  In  this 
case  one  can  consider  potentials  consisting  solely  of  multipole  “masses”  of 
high  order  for  which  the  interaction  decays  very  rapidly.  This  can  occur  if 
the  masses  (charges)  are  so  arranged  that  the  lower  order  moments  cancel. 
Thus,  the  effects  of  these  multipoles  can  be  ignored  beyond  a  relatively  short 
distance.  Wavelets  can  be  viewed  heuristically  as  high  order  multipoles. 

The  use  of  basis  functions  with  localized  support  provides  another  ben¬ 
efit:  the  matrices  for  estimation  of  the  model  parameters  will  be  sparse  (up 
to  any  finite  but  arbitrary  accuracy),  which  gives  us  a  chance  to  develop 
fast  algorithms  for  their  inversion.  The  coefficients  of  these  models  will  be 
estimated  directly  from  the  observed  data. 

3.1  Multiresolutior  QR  Algorithm 

Let  us  describe  our  approach  to  the  estimation  problem  in  more  detail.  As 
it  is  well  known,  the  least-squares  solution  to  a  system  of  linear  equations  is 
characterized  by  the  following 

Theorem:  Let  A  be  a  real  m  x  n  matrix  and  b  a  vector  of  length  m.  If  x 
satisfies 

AtAx  =  At  b 

then,  for  any  vector  y, 

||b  —  Ax||2  <  ||b-Ay||2. 

It  is  typical,  in  practical  estimation  problems,  for  the  matrix  product  AT A  to 
be  badly  conditioned.  If  this  is  the  case,  or  if  the  matrix  product  has  a  non- 
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trivial  null  space,  then  there  is  more  than  one  solution  to  the  least-squares 
problem.  Using  a  standard  QR  factorization,  it  is  typical  to  obtain  that  so¬ 
lution  to  the  least-squares  problem  that  minimizes  the  L2-norm.  However, 
for  the  real-world  estimation  problem,  such  solution  may  be  meaningless. 
We  introduce  a  new  approach,  based  on  a  physical  heuristic,  namely  that 
the  solution  we  desire  is  that  with  the  minimum  amount  of  detail.  That  is, 
we  attempt  to  minimize  the  high  frequency  content  of  the  solution.  Such 
procedure  is  feasible  in  a  multiresolution  setting,  where  frequency  content  is 
conveniently  split  across  several  scales,  allowing  us  to  work  with  each  fre¬ 
quency  range  separately.  The  Multiresolution  Rank  Revealing  QR  algorithm 
we  introduce  is  based  on  a  standard  Gramm-Schmidt  QR  factorization  but 
with  two  new  features:  (1)  unknowns  are  the  wavelet  coefficients,  and  (2) 
pivoting  is  restricted  to  appropriate  subspaces. 

As  an  example,  consider  the  following  system  of  linear  equations, 

Ax  =  y ,  (3) 

where  A  is  a  real  (m  x  n)  matrix,  and  we  explicitly  assume  m  <  n.  Let  us 
express  this  equation  in  the  following  equivalent  form, 

Ax  =  y ,  (4) 

where  A  =  AWT,  x  =  Wx,  and  W  denotes  the  wavelet  transform.  Since 
the  transform  is  orthogonal,  the  product  WTW  is  an  identity.  The  operation 
AWT  performs  wavelet  decomposition  along  rows  of  the  matrix  A,  and  thus 
the  columns  are  effectively  organized  into  subbands  corresponding  to  different 
scales.  We  next  construct  the  factorization 

AE  =  QR ,  (5) 

where  Q  is  an  orthogonal  ( m  x  r)  matrix,  R  is  an  upper  triangular  (r  x 
n)  matrix,  and  r  is  the  (numerical)  rank  of  A.  Note  that  r  <  m.  By 
construction,  we  have  QTQ  =  IT ,  which  denotes  the  (r  x  r )  identity  matrix. 
The  matrix  E  in  (5)  is  a  permutation  matrix,  and  serves  the  purpose  of 
re-ordering  the  columns  in  A  so  that  the  r  columns  in  Q,  which  form  an 
orthonormal  set,  are  formed  from  those  columns  in  A  which  are  largest  in 
/2-norm. 

The  solution  to  (4)  is  now  obtained  by  defining 

x  =  ERQt  y,  (6) 
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where  R  is  chosen  to  satisfy  RR  =  Ir.  Since  the  row  dimension  of  R  is  strictly 
less  than  the  column  dimension,  there  is  more  than  one  possible  choice  for 
R ,  and  we  describe  our  choice  below.  Substituting  (6)  into  (4),  we  have 

Ax  =  (QREt){ERQt)  y 
=  Q(RR)QTy 
=  ( QQT)y  ■ 

If  r  =  m,  then  QQT  is  an  identity,  and  we  have  y  -  Ax  =  0.  If  r  <  m,  then 
QQT  is  a  projection,  and  we  have 

||y-ix||  <  e 

if  y  lies  in,  or  nearly  in,  the  appropriate  subspace. 

Let  us  now  describe  the  procedure  for  constructing  Q.  To  begin  with,  we 
work  only  with  the  columns  of  averages  on  the  coarsest  scale  in  the  trans¬ 
formed  matrix  A.  From  among  these,  we  choose  the  column  with  largest 
/2-norm,  and  denote  this  column  vector  by  a^,  where  j\  is  the  actual  column 
index.  The  first  column  of  Q ,  denoted  qi,  is  obtained  as 

qi  IKll ' 

Thus,  qfqi  =  1,  and  we  now  multiply  all  remaining  columns  of  A  by  (I  — 
qiqfV  which  is  projection  onto  the  orthogonal  complement  of  {qi}.  The 
resub  cf  this  operation  is  that  all  remaining  columns  (excepting  a 2l)  in  the 
modified  matrix  A  are  orthogonal  to  qi. 

Now  repeat  this  step,  again  considering  only  columns  within  the  subband 
containing  averages  on  the  coarsest  scale,  and  choosing  from  among  these 
the  one  with  largest  /2-norm,  which  we  denote  by  aj2.  Obtain 


then  multiply  all  remaining  columns  of  A  by  (/-q2q^),  so  that  all  remaining 
columns  of  A  are  orthogonal  to  {qi,q2}-  We  repeat  this  process  until  one 
of  two  things  happens.  Either  we  have  “used  up”  all  vectors  in  the  subband 
of  averages,  or  the  remaining  vectors  in  this  subband  are  smaller  in  norm 
than  the  prescribed  cutoff  value.  In  either  case,  we  then  go  to  the  subband 
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containing  column  vectors  of  wavelet  coefficients  on  the  coarsest  scale,  and 
repeat  this  procedure  When  we  have  either  “used  up”  all  vectors  in  this 
subband,  or  the  vectors  remaining  in  this  subband  are  all  smaller  in  norm 
than  the  prescribed  cutoff,  then  we  move  to  the  subband  which  contains 
column  vectors  of  wavelet  coefficients  on  the  next  finer  scale.  This  procedure 
is  repeated  until  all  remaining  columns  in  A ,  which  are  being  overwritten  at 
each  step  of  the  process  with  their  projection  onto  the  orthogonal  complement 
of  the  column  vectors  in  Q,  are  smaller  in  norm  than  the  prescribed  cutoff 
value.  At  this  point,  we  have  r  vectors  in  Q  which  form  an  orthonormal  basis 
(to  within  prescribed  precision)  for  the  column  space  of  A. 

The  matrix  R  contains  coefficients  which  represent  the  columns  of  the 
product  AE  with  respect  to  the  orthonormal  basis  in  Q.  In  particular,  the 
jth  column  of  R  contains  the  coefficients  that  represent  the  jth  column  of 
AE  with  respect  to  the  orthonormal  basis  {qi, . . .  ,qr},  and  by  construction 
R  is  upper  triangular.  To  construct  R,  so  that  RR  =  Jr,  we  invert  the  leading 
r  columns  of  R  which,  taken  together,  form  an  upper  triangular  (r  x  r)  block, 
and  thus  obtain  the  first  r  rows  of  R.  We  then  set  the  remaining  (n  -  r ) 
rows  of  R  to  zero.  This  has  the  effect  of  selecting  only  r  non-zero  entries 
in  the  solution  vector  x,  while  setting  the  remaining  entries  to  zero.  The 
algorithm  automatically  selects  which  entries  in  x  correspond  to  the  most 
heavily  weighted  columns  of  A,  which  for  physical  data  should  correspond 
to  the  lower  frequency  features  of  the  right-hand  side.  We  also  emphasize 
the  selection  of  lower  frequency  components  by  beginning  the  procedure  with 
columns  on  the  low  frequency  side  and  proceeding  from  lower  to  higher.  . 

Figures  1-3  on  the  following  pages  compare  results  of  Multiresolution 
QR  and  a  standard  QR  for  the  following  interpolation  problem.  In  (3),  the 
right-hand  vector  y  represents  values  of  a  function,  say  y(t),  at  a  set  of 
points  {fi, . . .  ,*m}.  Thus,  y  =  {yi, . . .  ,ym},  where  y,-  =  y{U).  We  build  an 
expansion  using  multiwavelet  scaling  functions  on  some  scale,  which  takes 
the  form 

n— 1 

y(t)  =  J2  xj+Mnt  -  j)  ■ 

j= 0 

Imposing  the  condition  Y(U)  =  y*,  leads  to  the  system  of  equations 

n— 1 

-  j)  =  yi} 
j= o 

and  comparing  with  (3)  we  see  that  the  entry  of  the  matrix  A  is  given 
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by  (f>(nti  -  j),  and  the  unknowns  {xi,...,x„}  are  the  coefficients  of  the 
expansion. 

We  choose  the  scale  of  the  expansion  so  that  the  number  n  of  coefficients 
on  the  finest  scale  is  much  larger  than  the  number  m  of  data  points.  How¬ 
ever,  the  algorithm  will  select  a  smaller  number  of  coefficients  to  represent 
the  final  solution.  If  the  data  is  smooth,  then  the  matrix  A  will  have  low 
numerical  rank  at  the  required  precision,  and  the  final  representation  will  be 
correspondingly  sparse. 

The  next  three  pages  contain  examples  which  compare  the  results  of  the 
approached  outlined  above  with  results  of  the  usual  least  squares  solution 
method. 
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Figure  3.  The  usual  QR  solution  (top)  and  multiresolution  QR  (bottom) 
for  the  least-squares  problem  with  a  null  space. 
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4  The  Next  Generation  of  Gravity  Models 

The  choice  of  basis  function  used  to  represent  the  gravity  field  is  of  enormous 
importance  since  this  choice  has  the  single  greatest  influence  on  the  structure 
and  performance  of  algorithms  that  will  be  used  for  evaluation  and  estimation 
of  the  model.  Under  this  grant  we  have  developed  a  new  type  of  basis  for 
estimation  of  gravity  fields,  though  we  have  not  yet  developed  the  actual 
models,  and  we  explain  below  the  main  points  of  the  new  construction.  We 
foresee  that  the  next  generation  of  gravity  models  will  use  such  bases  in 
conjunction  with  the  algorithms  discussed  above. 

4.1  Optimal  Bases  on  an  Interval 

Let  us  first  provide  a  brief  description  of  the  prolate  spheroidal  wave  functions 
(PSWFsj  introduced  by  Slepian  et.  al.  [20,  16]. 

The  PSWFs  are  defined  as  the  eigenfunctions  of  the  operator  Fc,  where 
Fc  is  defined  by 

Fc{<t>){x)  =  (7) 

where  Fc  :  L 2  [-1, 1]  ->■  L2  [-1, 1],  and  c  is  a  positive  real  constant  (band- 
limit).  The  PSWFs  satisfy 

A Mx)  =  j\icxtiPj(t)dt,  (8) 

where  the  eigenvalues  A j,  j  =  0, 1, 2, . . .,  are  all  non  zero  and  simple,  and  are 
arranged  so  that  |Aj|  >  |AJ+1|.  The  eigenvalues  A  j  are  either  real  or  pure 
imaginary  depending  on  the  parity  of  the  eigenfunction  ipj.  These  eigenfunc¬ 
tions  are  also  eigenfunctions  of  the  operator  Qc  =  ^F*FC  and  satisfy 

=  (9) 

with  eigenvalues 

^  =  ^=0’1'2 .  (10) 

For  large  c  the  first  approximately  2c/ tt  eigenvalues  Hj  are  close  to  1.  These 
are  followed  by  O(logc)  eigenvalues  which  decay  exponentially  fast  from  1  to 
almost  zero.  The  remaining  eigenvalues  are  very  close  to  zero. 
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In  addition,  there  exists  a  strictly  increasing  sequence  of  real  numbers 
Vo  <  Vi  <  %  <  -  • '  such  that  the  functions  ipj  in  (8)  are  eigenfunctions  of  the 
following  differential  operator  [20], 

Lipj  =  ^-(1  -  x2)  +  2x +  c2x2^j  ipj(x)  =  VjTpj(x)  •  (11) 

The  eigenfunctions  of  L  have  been  known  as  the  angular  prolate  spheroidal 
functions  before  the  connection  with  (8)  was  discovered  [20].  We  note  that, 
in  the  limit  c  — >  0,  it  follows  from  (11)  that  ipj  become  the  Legendre  poly¬ 
nomials. 

For  any  n  >  0,  the  first  n  functions  ipj,  j  =  0, . . . ,  n— 1,  form  a  Chebyshev 
system  [14,  15].  In  particular,  the  number  of  zeros  of  ipj  in  [—1, 1]  is  equal 
to  j. 

Although  the  functions  ipj  are  defined  on  the  interval,  they  are  easily 
extended  to  the  whole  line  using  the  right  hand  side  of  (9)  as  the  definition 
of  the  extension.  The  functions  ipj  are  orthogonal  on  both  the  interval  [—1,1] 
and  the  real  line  (—00,00),  and  we  set 

J  ^  ipj(x)  ipi{x)  dx  =  5ji ,  (12) 

and 

/OO  j 

ipj(x)  ipi(x)  dx  =  —8ji .  (13) 

-OO  /ij 

We  note  that  in  the  original  papers  [20,  16,  19]  the  functions  are  chosen  to 
be  orthonormal  on  (—00,00). 

The  definition  (8)  implies  that 

eicxt  _  A jipj(x)ipj(t) ,  (14) 

j= 0 

and  if  we  keep  approximately  2c/tt  +  K  log  c  terms  in  (14),  where  K  —  K(e) 
is  a  constant,  we  obtain  a  close  approximation  to  elcxt  for  any  positive  e.  This 
is  the  most  economical  expansion  of  this  type  for  the  exponential. 

The  PSWFs  have  been  used  in  signal  processing  for  some  time,  especially 
the  first  function,  ip0(x ),  since  it  provides  the  optimal  window  for  a  given 
bandwidth  in  terms  of  concentration  in  the  time-frequency  domain.  Yet, 
their  use  has  not  been  wide.  In  the  next  section  we  describe  several  new 
developments  that  will  provide  a  path  for  a  wider  use  of  these  functions  in 
signal  processing  and  numerical  analysis. 
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4.2  Generalized  Gaussian  Quadratures  for  Exponen¬ 
tials 

The  generalized  Gaussian  quadratures  for  exponentials  has  been  developed 
recently  [22,  6].  Within  the  tirst  approach  [22],  the  authors  construct  the 
generalized  Gaussian  quadratures  for  the  prolate  spheroidal  wave  functions 
using  the  fact  that,  for  any  n,  the  first  n  of  these  functions  form  a  Chebyshev 
system  [14,  15].  For  a  given  accuracy  t  and  a  choice  of  n,  if  follows  from 

(14)  that  such  quadratures  are  also  valid  for  exponentials.  Alternatively,  a 
new  type  of  the  generalized  Gaussian  quadratures  for  exponentials  has  been 
obtained  directly  [6].  These  quadratures  are  parameterized  by  eigenvalues 
of  a  Toeplitz  matrix  which  is  constructed  from  the  trigonometric  moments 
of  a  positive  measure.  For  a  given  accuracy  e,  selecting  an  eigenvalue  close 
to  t  yields  an  approximate  quadrature  for  that  accuracy.  These  quadratures 
can  be  used  to  approximate  and  integrate  other  essentially  bandlimited  func¬ 
tions  such  as,  for  example,  Bessel  functions  or  the  prolate  spheroidal  wave 
functions. 

Let  us  define  the  bandlimited  functions  with  the  bandlimit  c  as  a  class  of 
functions  that  can  be  represented  via  a  linear  combination  of  exponentials  of 
the  form  exp  (i bx)  with,  for  example,  ^-bounded  coefficients,  where  b  is  any 
real  number  such  that  \b\  <  c. 

It  turns  out  that,  for  any  accuracy  e  >  0  and  any  bandlimit  c  >  0,  there  is 
a  set  of  M  functions,  {exp  (ictkx)}%Lv  where  the  nodes  tk  =  tk(e,  c),  |t*|  <  1 
and  the  coefficients  ak  =  otk(b)  are  such  that 

M 

|  exp  (i bx)  -  ^  ak(b)  exp  (icifcar)|  <  e .  (15) 

k=l 

The  set  of  functions  {exp  (\ctkx)}^=l  can  be  viewed  as  an  approximate  basis. 

In  order  to  find  the  nodes  {tk}  in  (15),  we  solve  the  following  problem 
[6],  described  here  in  a  slightly  more  general  setting  than  is  needed  to  obtain 

(15) .  Let  us  consider  integrals  of  the  form 

u(x)  =  i:  exp  (2  ctx)dfi(t),  (16) 

where  dfj,(t)  is  a  measure.  Typically,  dpi(t)  =  w(t )  dt,  where  w  is  a  weight 
function,  that  is  w  is  a  real,  non-negative,  integrable  function  with  /ix  w(r)dr  > 
0.  To  obtain  (15)  the  weight  is  defined  as  w  =  1. 
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For  a  given  bandlimit  2c  >  0  and  accuracy  e2  >  0,  we  approximate  u(x) 
on  the  interval  [—1,1]  using  the  sum 

M 

u(x)  exP  ( 2ctk  *),  (17) 

jt=i 

where  Wk  >  0  and  M  =  M(c,  e2),  so  that 

|u(x)  —  u(:r)|  <  e2  for  x  €  [—1, 1] .  (18) 

The  number  of  terms,  M,  is  optimal.  Solving  this  problem  involves  finding 
the  eigenvalues  and  eigenvectors  of  the  Toeplitz  matrix  constructed  using  the 
values  of  u(x )  discretized  at  the  equally  spaced  nodes  and  interpreted  as  the 
trigonometric  moments  of  a  positive  measure  [6]. 

Once  the  nodes  are  computed,  the  set  of  functions  {exp  (ict^a;)}^!  can 
serve  as  an  approximate  basis  on  the  interval  [—1,1]  in  (15).  Such  bases  can 
be  organized  into  a  hierarchical  structure  similar  to  multiwavelets.  We  will 
report  these  results  elsewhere. 

The  representation  in  (15)  retains  the  property  of  disjoint  support  similar 
to  that  of  a  multiwavelet  basis.  On  the  other  hand,  it  requires  significantly 
fewer  terms  than  the  representation  with  orthogonal  polynomials.  Also,  one 
can  think  of  the  bandlimit  c  as  an  analogue  of  the  degree  in  the  case  of 
polynomials  and,  unlike  in  that  case,  there  is  no  constraint  on  the  bandwidth 
c.  This  is  because  the  distance  between  the  nodes  is  0(l/c)  and,  thus, 
the  quadrature  nodes  do  not  significantly  concentrate  near  the  ends  of  the 
interval  as  do,  for  example,  the  Legendre  nodes.  These  properties  of  the 
representations  using  exponentials  lead  to  a  number  of  new  algorithms  that 
are  being  developed  and  will  be  presented  elsewhere. 


5  Spectral  Deferred  Corrections  for  Solution 
of  Ordinary  Differential  Equations 

This  is  a  relatively  new  method  for  solution  of  ordinary  differential  equations 
[12],  which  shows  great  promise  for  further  optimizing  the  computation  of 
satellite  ephemerides.  The  basic  idea  is  to  first  obtain  a  crude  approximation 
using  a  simple  first-order  method,  for  example  an  explicit  Euler  scheme,  then 
to  iteratively  reduce  (correct)  the  error.  Methods  based  on  this  approach 
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are  known  collectively  as  deferred  correction  methods,  and  have  been  in 
use  for  some  time.  However,  in  [12]  some  improvement  over  the  classical 
methods  is  introduced,  namely  the  use  of  orthogonal  polynomials  to  represent 
the  solution,  together  with  spectral  integration,  which  greatly  increases  the 
stability  of  high-order  deferred  correction  methods.  This  approach  leads 
to  methods  of  almost  arbitrarily  high  order  which  require  a  relatively  low 
number  of  function  calls  in  their  implementation. 

We  have  constructed  a  preliminary  version  of  a  spectral  deferred  correc¬ 
tion  code  for  the  orbit  propagation  problem,  and  initial  tests  show  promise. 
Work  is  in  progress  towards  developing  a  final  version  which  is  capable  of 
replacing  current  operational  orbit  propagation  codes. 

6  Transfer  of  Technology 

The  latest  version  of  the  B-spline  cube  model  has  been  transferred  to  Space 
Warfare  Command  (SWC)  in  Colorado  Springs,  to  our  point (s)  of  contact, 
Joseph  Liu  and  Mark  Storz,  both  of  SWC.  A  declassified  copy  of  currently 
operational  code  for  solving  ODES  has  been  transferred  to  us,  which  allows 
us  to  make  direct  comparisons  with  new  codes  we  are  considering. 


References 

[1]  B  Alpert.  A  class  of  bases  in  L 2  for  the  sparse  representation  of  integral 
operators.  SIAM  J.  Math.  Anal,  24(l):246-262,  1993. 

[2]  B.  Alpert,  G.  Beylkin,  R.R.  Coifman,  and  V.  Rokhlin.  Wavelet-like  bases 
for  the  fast  solution  of  second-kind  integral  equations.  SIAM  J.  Sci. 
Statist.  Comput.,  14(1):159-174,  1993.  Technical  report,  Department  of 
Computer  Science,  Yale  University,  New  Haven,  CT,  1990. 

[3]  B.  Alpert,  G.  Beylkin,  D.  Gines,  and  L.  Vozovoi.  Adaptive  solution  of 
partial  differential  equations  in  multiwavelet  bases.  Technical  Report 
APPM  409,  Univ.  of  Colorado  at  Boulder,  1999. 

[4]  G.  Beylkin.  On  the  representation  of  operators  in  bases  of  compactly 
supported  wavelets.  SIAM  J.  Numer.  Anal.,  29(6) :  1716— 1740,  1992. 


22 


[5]  G.  Beylkin.  On  wavelet-based  algorithms  for  solving  differential  equa¬ 
tions.  In  John  J.  Benedetto  and  Michael  W.  Frazier,  editors,  Wavelets: 
Mathematics  and  Applications,  pages  449-466.  CRC  Press,  1994. 

[6]  G.  Beylkin  and  L.  Monzon.  On  generalized  Gaussian  quadratures  for 
exponentials  and  their  applications.  Technical  Report  APPM  452,  Univ. 
of  Colorado  at  Boulder,  2000.  Submitted  for  publication. 

[7]  A.  Cohen,  I.  Daubechies,  B.  Jawerth,  and  P.  Vial.  Multiresolution  anal¬ 
ysis,  wavelets  and  fast  algorithms  on  an  interval.  Comptes  Rendus  Acad. 
Sc.  Paris,  1992. 

[8]  A.  Cohen,  I.  Daubechies,  and  P.  Vial.  Wavelets  on  the  interval  and 
fast  wavelet  transforms.  Applied  and  Computational  Harmonic  Analysis, 
1(1):54-81,  1993. 

[9]  R.R.  Coifman  and  Y.  Meyer.  Remarques  sur  l’analyse  de  Fourier  a 
fenetre.  C.  R.  Academie  des  Sciences,  312(1):259 — 261,  1991. 

[10]  I.  Daubechies.  Orthonormal  bases  of  compactly  supported  wavelets. 
Comm.  Pure  Appl.  Math.,  41(7):909-996,  1988. 

[11]  Department  of  Defense  World  Geodetic  System  1984.  Defense  Mapping 
Agency  Technical  Report,  DMA  TR  8350.2,  1987. 

[12]  A.  Dutt,  L.  Greengard,  and  V.  Rokhlin.  Spectral  deferred  correction 
methods  for  ordinary  differential  equations.  BIT 40(2):241-266,  2000. 

[13]  P.  Federbush.  A  mass  zero  cluster  expansion.  Communication  Math. 
Phys.,  81:327-340,  1981. 

[14]  S.  Karlin  and  W.  J.  Studden.  Tchebyche ff  systems:  With  applications 
in  analysis  and  statistics.  Interscience  Publishers  John  Wiley  &  Sons, 
New  York-London-Sydney,  1966.  Pure  and  Applied  Mathematics,  Vol. 
XV. 

[15]  M.  G.  KreTn  and  A.  A.  Nudel’man.  The  Markov  moment  problem  and 
extremal  problems.  American  Mathematical  Society,  Providence,  R.I., 
1977.  Ideas  and  problems  of  P.  L.  Cebysev  and  A.  A.  Markov  and  their 
further  development,  Translations  of  Mathematical  Monographs,  Vol. 
50. 


23 


[16]  H.  J.  Landau  and  H.  O.  Poliak.  Prolate  spheroidal  wave  functions, 
Fourier  analysis  and  uncertainty  II.  Bell  System  Tech.  J.,  40:65-84, 
1961. 

[17]  P.-G.  Lemarie  and  Y.  Meyer.  Ondelettes  et  bases  hilbertiennes.  Revista 
Mat.  Iberoamericana,  2:1-18,  1986. 

[18]  H.  S.  Malvar.  Lapped  Transforms  for  Efficient  Transform/Subband  Cod¬ 
ing.  IEEE  Trans.  Acoust.,  Speech,  Signal  Processing ,  38(6):969-978, 
1990. 

[19]  D.  Slepian.  Prolate  spheroidal  wave  functions,  Fourier  analysis  and 
uncertainty  IV.  Extensions  to  many  dimensions;  generalized  prolate 
spheroidal  functions.  Bell  System  Tech.  J .,  43:3009-3057,  1964. 

[20]  D.  Slepian  and  H.  O.  Poliak.  Prolate  spheroidal  wave  functions,  Fourier 
analysis  and  uncertainty  I.  Bell  System  Tech.  J.,  40:43-63,  1961. 

[21]  J.  O.  Stromberg.  A  modified  franklin  system  and  higher-order  spline 
systems  on  Rn  as  unconditional  bases  for  Hardy  spaces.  In  Conference 
in  harmonic  analysis  in  honor  of  Antoni  Zygmund,  Wadworth  math, 
series,  pages  475-493,  1983. 

[22]  H.  Xiao,  V.  Rokhlin,  and  N.  Yarvin.  Prolate  spheroidal  wave  functions, 
quadrature,  and  interpolation.  Technical  Report  YALEU/DCS/RR- 
1199,  Yale  University,  2000. 


24 


